Comprehensive analysis of m6A modification lncRNAs in high glucose and TNF-α induced human umbilical vein endothelial cells

N6-methyladenosine (m6A) RNA methylation, as a reversible epigenetic modification of mammalian mRNA, holds a critical role in multiple biological processes. m6A modification in Long non-coding RNAs (lncRNAs) has increasingly attracted more attention in recent years, especially in diabetics, with or without metabolic syndrome. We investigated via m6A-sequencing and RNA-sequencing the differentially expressed m6A modification lncRNAs by high glucose and TNF-α induced endothelial cell dysfunction in human umbilical vein endothelial cells. Additionally, gene ontology and kyoto encyclopedia of genes and genomes analyses were performed to analyze the biological functions and pathways for the target of mRNAs. Lastly, a competing endogenous RNA network was established to further reveal a regulatory relationship between lncRNAs, miRNAs and mRNAs. A total of 754 differentially m6A-methylated lncRNAs were identified, including 168 up-regulated lncRNAs and 266 down-regulated lncRNAs. Then, 119 significantly different lncRNAs were screened out, of which 60 hypermethylated lncRNAs and 59 hypomethylated lncRNAs. Moreover, 122 differentially expressed lncRNAs were filtered, containing 14 up-regulated mRNAs and 18 down-regulated lncRNAs. Gene ontology and kyoto encyclopedia of genes and genomes analyses analyses revealed these targets were mainly associated with metabolic process, HIF-1 signaling pathway, and other biological processes. The competing endogenous RNA network revealed the regulatory relationship between lncRNAs, miRNAs and mRNAs, providing potential targets for the treatment and prevention of diabetic endothelial cell dysfunction. This comprehensive analysis for lncRNAs m6A modification in high glucose and TNF-α-induced human umbilical vein endothelial cells not only demonstrated the understanding of characteristics of endothelial cell dysfunction, but also provided the new targets for the clinical treatment of diabetes. Private information from individuals will not be published. This systematic review also does not involve endangering participant rights. Ethical approval will not be required. The results may be published in a peer-reviewed journal or disseminated at relevant conferences.


Introduction
Type 2 diabetes mellitus (T2DM), is a chronic endocrine and metabolic disorder syndrome, accounting for more than 90% of diabetic patients. [1,2] It is characterized by hyperglycemia, hyperlipidemia, hypertension and insulin resistance. T2DM, which is common along with complications such as diabetic nephropathy, retinopathy, and vascular disease. [3,4] Accumulating evidence suggested that patients with T2DM have an increased risk of vascular complications. Meanwhile, the development of coronary heart disease, cardiomyopathy, arrhythmias and peripheral artery disease are the primary cause of mortality in T2DM. [5,6] Owing to its complex and diverse pathogenesis, the mechanism of the impairment of cardiovascular disease in diabetes has not been fully established. However, oxidative stress mediated by hyperglycemia and activation of inflammatory damage have been recognized as key underlying events. In addition, many clinical studies have demonstrated that vascular endothelial cell proliferation, basement membrane thickening, and deposition of hyaline-like substances may contribute to the mechanism of diabetic vascular disease. [7,8] LS and MG contributed equally to this article. Long non-coding RNAs (lncRNAs) are defined as RNA molecules with a transcript length of more than 200 nt. Several studies have shown that lncRNAs participate in many life activities, such as dosage compensation effect, epigenetic regulation, cell cycle and differentiation regulation. [9,10] Abnormal expression of lncRNAs has been found to play an important role in multiple diseases, especially in T2DM and its complications. Although it has gained increasingly deep research, the underlying molecular mechanism of lncRNAs was still mostly uncharacterized in diabetic angiopathies.
N6-methyladenosine (m6A), the most prevalent and abundant internal modification of RNA in eukaryotic cells, which plays essential roles in mRNA metabolism and multiple biological processes. [11,12] The dysregulation of m6A modification is associated with almost every stage of RNA metabolism, ranging from RNA splicing, nuclear export and translation to stability. [13] Current studies have identified that aberrant m6A methylation in diabetic angiopathies, including heart failure, vascular calcification and pulmonary hypertension. [14][15][16] However, the molecular mechanism of m6A methylation in multifarious diseases is still not fully understood. Therefore, proclaiming the potential of m6A machinery as novel targets for prevention and treatment diabetic angiopathies will attract more attention.
In addition, accumulating research has demonstrated that abundant lncRNAs are also highly modified with m6A to perform their functions. Whereas the pathogenic mechanism of m6A modification lncRNAs in diabetic atherosclerosis have not been comprehensively clarified. In addition, TNF-α as proinflammatory cytokines is an important trigger of endothelial cell dysfunction. Meanwhile, it has been demonstrated they play a crucial role in promoting endothelial cell dysfunction via diverse mechanisms, such as oxidative stress and inflammatory response. Accordingly, high glucose and TNF-α induced endothelial cell dysfunction model in human umbilical vein endothelial cells (HUVECs) was established to reveal the function of m6A modification lncRNAs in atherosclerosis. It also provides a reference for exploring new targets for the diagnosis and treatment of atherosclerosis.

Cell culture and pretreatment
HUVECs were purchased from iCell Biological Technology (Shanghai, China). HUVECs cells were cultured in Dulbecco's modified Eagle's medium (DMEM) supplemented with 10% fetal bovine serum (FBS; BI), 100 U/mL penicillin, 100 mg/mL streptomycin at 37°C in 5% CO 2 -humidified incubator. To estimate the effects of high glucose on HUVECs, the cells were pretreated with medium containing 5.5 mM glucose, and then subjected to high glucose and TNF-α treatment. Meanwhile, 25 mM mannitol was used as an osmolarity control condition. Lastly, culture medium concentration was generated by adding 25 mM glucose and 5 ng/mL TNF-α as the final qualification. [17] Afterwards, different treatment cells in two groups were cultured in an incubator for 48 hours for the further experiment.

MeRIP and RNA library preparation and sequencing
Total RNA was extracted using TRIzol reagent (Invitrogen, CA), according to the manufacturer's protocol. RNA integrity and quality were analyzed via nanodrop and gel electrophoresis. m6A RNA immunoprecipitate was performed by the GenSeqTM m6A-MeRIP Kit (GenSeq Inc., Malaysia), according to the manufacturer's instructions. Both the input samples without immunoprecipitation and the m6A samples with immunoprecipitation were used for RNA-seq library generation. The library quality was estimated with Bioptic Qsep100 Analyzer (Bioptic lnc., Taiwan, China). MeRIP-Seq service was rendered by Shanghai Biotechnology Corporation (Shanghai, China). Library sequencing was transacted on an illumina NovaSeq 6000 with 150bp paired-end reads.

Data analysis of MeRIP and RNA sequencing
The quality of the original sequencing data was evaluated through FastQC software (v0.11.7). The raw reads were trimmed using Cutadapt (v2.5) and HISAT2 software (v2.1.0) was aligned to the Ensembl database (GRCh38/ hg38). Cutadapt (v2.5) was used to trim adapters and filter for sequences, remaining reads were aligned to the human Ensemble database (GRCh38/hg38). The MeRIP enriched peaks were identified using exomePeak (v2.13.2). Differentially m6A-methylated lncRNAs peaks between the control and model group were analyzed using exomepeak software by Poisson-Gamma performed. Identified m6A peaks that P < .05 were chosen for the denovo motif analysis using HOMER (v4.10.4).
The gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) enrichment analyses were performed for the differentially methylated related genes. The competing endogenous RNA network was constructed through the Cytoscape software (v3.6.1). The targets of lncRNAs were predicted from the LnaACTdb and LncTarD tools.

Characteristics of m6A methylation of lncRNAs
A total of 7581 m6A peaks were detected within 557 lncRNAs in control group, while 5983 m6A peaks were recognized within 459 lncRNAs in model group ( Fig. 1A and B). In addition, 5461 m6A peaks and 65 lncRNAs were identified in both groups. Incidentally, 4937 m6A peaks were all shared in the two groups, including 2492 up-regulated peaks and 2445 down-regulated peaks (Fig. 1C). Meanwhile, 434 differently m6A methylated lncRNAs were screened, containing 168 up-regulated lncRNAs and 266 down-regulated lncRNAs (Fig. 1D). Overall, these results indicated that the degree of the m6A modification in lncRNAs was higher in the control group than in the model group.
Inasmuch to further expose the positional relationship of m6A methylated lncRNAs, we divided into the following categories, containing coding sequences, 3ʹ-untranslated regions (UTRs), 5ʹ-UTRs, and exon. In control group, m6A methylated lncRNAs levels were increased in the exon (93.350% vs 92.25%) and 3ʹ-UTRs (3.80% vs 3.21%), when compared to the model group (Fig. 2C). On the contrary, m6A methylated lncRNAs levels were decreased in coding sequences when compared to the model (2.85 % vs 4.28 %). Moreover, approximately 0.27% of the m6A methylated lncRNAs were only appeared in the 5ʹ-UTRs region of the model group (Fig. 2D).

Abundance of m6A peaks and conserved m6A motifs in lncRNAs
Regarding the abundance of the m6A peaks in lncRNAs, we found that 77.13% of the lncRNAs in the control group contained m6A peaks, which appeared marginally more than the unimodal value calculated at 75.86% in the model group. The respective percentages comparing different numbers of peaks were also determined with two peaks, three peaks, and more than three peaks being 15.81 versus 16.66, 3.92% versus 5.10 and 3.14% versus 2.38%, respectively, for the control versus model group (Fig. 3A).
To analyze the conserved motif of m6A methylated lncRNAs, we selected the sequences of the peaks with the highest enrichment factor in two groups. The motif sequence was compared with the peaks with the highest enrichment ratio of lncRNA. It was found that GGACCG sequence was one of the conserved motif sequences of lncRNA based on P value (Fig. 3B).

Synopsis of differentially expressed lncRNAs in HUVECs
As shown in Figure 4A, the results indicated that these lncRNAs have different expression patterns in the two groups. There were 32 significant differentially expressed lncRNAs in HUVECs, which were 14 up-regulated and 266 down regulated (Fig. 4B). The top 10 up-regulated and down-regulated lncRNAs are listed in Table 2.

Association of m6A methylation and expression of lncRNAs
Combining the results of methylation sequencing and RNA sequencing, we found that 185 lncRNAs were hypermethylated, of which was only 4 differentially expressed lncRNAs were up-expressed. A total of 292 lncRNAs were hypomethylated, among them, 13 lncRNAs were down-expressed (Fig. 5A).
To analyze the correlation between lncRNAs methylation level and expression level, a correlation graph was constructed using the fold enrichment of lncRNA m6A methylation and expression value in terms of FPKM. The results demonstrated that there was a statistically significant positive correlation between methylation and expression levels of lncRNAs in the control and model groups (Fig. 5B).

Functional of differentially m6A methylation and expression of lncRNAs
Owing to manifest the role of differentially methylated and expression of lncRNAs in the occurrence and development of endothelial cell dysfunction, GO and KEGG pathway analyses were performed on the genes located near differentially lncRNAs. GO results revealed that these genes primarily participate in the metabolic process, positive regulation of biological process, biological regulation, and cellular process in the biological process category. In terms of cellular component, genes are associated with cell part, organelle part, and so on. In terms of molecular function, genes are primarily involved in contains binding (Fig. 6A). Meanwhile, we have demonstrated that five lncRNAs MEG3, MALAT1, FTX, XIST, and NEAT1  were mainly concentrated on the GO analysis, which deserves more attention in future research.
KEGG analysis found that the most important signaling pathways associated with genes were significantly enriched in PI3K-Akt signaling pathway, FOXO signaling pathway, pancreatic cancer and HIF-1 signaling pathway (Fig. 6B). These genes, such as PI3K, FOXO, HIF-1α, TGF-1β, MAPK1, EGFR, and BCL2 participated in the function of these pathways and accelerated the degree of endothelial cell dysfunction.

Construction of lncRNA-miRNA-mRNA network
Five significant lncRNAs were screened out from 51 differentially methylated and expression lncRNAs, of which were associated with endothelial cell dysfunction. The network of lncRNA-miRNA-mRNA was constructed by Cytoscape software. It consisted of the screened lncRNAs and combined with predicted corresponding miRNAs and mRNAs from the LnaACTdb and LncTarD software. The network contained 5 lncRNAs, 231 miRNAs, and 273 mRNAs (Fig. 7). It is evident from the competing endogenous RNA network that lncRNAs regulate the complex relationship between miRNAs and mRNAs. The network will contribute to understanding of the functions and mechanisms of lncRNAs in endothelial cell dysfunction.

Discussion
In this research, we have investigated the m6A modification lncRNAs by high glucose and TNF-α influenced in HUVECs. The results performed significantly different between control and model groups with the abundance and distribution  of m6A modification. In addition, the ratio of m6A modified lncRNAs in control group was remarkably greater than the model group, which indicated m6A modification down regulation of lncRNA expression in the model group. Meanwhile, previous research demonstrated that METTL14-mediated m6A methylation modification suppressed pyroptosis and diabetic cardiomyopathy by down-regulating TINCR. [18] Interestingly, the m6A methylation of XIST was prominently reduced, but the level of XIST was obviously over expression. [19] Notwithstanding, no studies have authenticated the absolute regulatory relationship between m6A methylation levels and lncRNA expression, which requires further in-depth research.
The targets of m6A methylation lncRNAs were analyzed through GO enrichment, including biological processes, molecular function and cellular components. Biological processes results were mainly involved in metabolic process, positive regulation of biological process, biological regulation, and cellular process. Molecular function enrichment analysis was mostly concentrated on binding. Meanwhile, cellular components consist of macromolecular complex, cell part, and organelle part. These biological processes were mostly  related to glucose and lipid metabolism, indicating that m6A modified lncRNAs perhaps involved in regulating insulin secretion and reducing blood glucose. GO analysis results revealed that majority of the m6A modified lncRNAs were hypomethylated sites, of which specific mechanisms still need further investigation.
KEGG pathways demonstrated that most of m6A modified lncRNAs were down-regulated in the model group. As an upstream pathway of FOXO, PI3K/Akt could reduce Akt phosphorylation and activate insulin secretion to reduce the level of blood glucose. Additionally, research suggests that regulating PI3K/Akt pathway alleviates oxidative stress and inhibits apoptosis in diabetic cardiomyopathy, together with attenuated diabetic vasculopathy. [20,21] FOXO signaling pathway could regulate insulin signaling, gluconeogenesis and immune cell migration in diabetes. Previous research has declared that metformin might alleviated HUVECs apoptosis and vascular endothelial injury by regulating FOXO protein. [22] Meanwhile, activated PI3K/Akt might regulate glycolysis through the HIF-1α pathway and reduce vascular damage under hypoxic conditions. [23] Altogether, these m6A modified lncRNAs mediated pathways could harmonize blood glucose, diminish the degree of atherosclerosis and alleviate endothelial cell dysfunction.
In addition, our results manifested that lncRNAs NEAT1, XIST, MALAT1, FTX, and MEG3 were differentially expressed in two groups, which might play a crucial role in m6A modification and endothelial cell dysfunction of HUVECs. NEAT1 has been reported to participate in multiple diabetic metabolic syndrome, which accelerates the occurrence and development of diabetic nephropathy by sponging miR-23c. [24] Others have claimed that lncRNA NEAT1 regulated diabetic retinal epithelial-mesenchymal transition via regulating miR-204/ SOX4 axis. [25] Our research demonstrated that the expression level of NEAT1 was significantly higher in the model group, which might be related to m6A modification. LncRNA XIST was downregulated in high glucose treated podocytes, accompanied with increased apoptosis of podocytes. Furthermore, lncRNA XIST protects podocyte from high glucose-induced cell injury by sponging miR-30 and regulating AVEN expression in diabetic nephropathy. [26] Additionally, m6A modification of lncRNAs levels has been observed in various prime aggravators of diabetic nephropathy pathogenesis. [27] The targets of XIST were involved in BCL2, STAT3, MAPK1, and VEGF, which mainly concentrate on the immune regulation, oxidative stress, inflammatory response.
What's more, current studies suggest that the lncRNAs NEAT1, XIST, and MALAT1 participate in the pathogenesis of T2DM and highlight their potential as diagnostic biomarkers, especially MALAT1, which is highly expressed in the serum of patients with coronary atherosclerosis heart disease, and it has high value in the diagnosis and prediction of in-stent restenosis. [28] Inhibition of MALAT1 has the potential to protect the retina from oxidative damage and to prevent diabetic retinopathy. [29] FTX has been reported to be associated with several tumor progressions, such as hepatocellular carcinoma, [30] renal cell carcinoma, [31] and colorectal cancer. [32] Nevertheless, we have revealed that it was associated with endothelial cell dysfunction, and further experiments are required and validate their effects and mechanisms. LncRNA MEG3 is related with multiple biological processes, containing proliferation, apoptosis and inflammation response. Previous research has illustrated MEG3 could alleviate high glucose inducing apoptosis and inflammation via inhibiting NF-κB pathway by targeting miR-34a/SIRT1 axis. [33] In addition, ROCK2, TNF, NOTCH1, and HIF-1α have also played crucial regulatory roles during atherogenesis. Accordingly, controlling blood glucose and reducing the stimulation of inflammatory factors have obvious promoting effects on reducing vascular damage and improving atherosclerosis. However, our study had a few limitations regarding the analysis of m6A-sequencing and RNA-sequencing. All the results were only based on association studies and bioinformatic analyses, which need further experiments to verify the results. In addition, the functions of m6A modification lncRNAs and their targets related to endothelial dysfunction require further investigation.

Conclusion
In summary, we firstly established a comprehensive analysis to investigate the m6A modification of lncRNAs in HUVECs to predict the mechanism of high glucose and TNF-α induced endothelial cell dysfunction. Despite the above limitations, our research still provides a new theoretical basis for the study of diabetic endothelial cell dysfunction, as well as new targets and reference for clinical treatment of atherosclerosis.